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ABSTRACT 

The Kepler mission has discovered more than 4000 exoplanet candidates. Many of t hem are in 
systems with tightly packed inner planets (STIPs). Inside-Out Planet Formation (lOPF) (Chatterjee 
& Tan||2014) has been proposed as a scenario to explain these systems. It involves sequeritial in situ 
planet formation at the local pressure maximum of a retreating dead zone inner boundary (DZIB). 
Pebbles accumulate at this pressure trap, which builds up a pebble ring, and then a planet. The 
planet is expected to grow in mass until it opens a gap, which helps to both truncate pebble accretion 
and also induce DZIB retreat that sets the location of formation of the next planet. This simple 
scenario may be modified if the planet undergoes significant migration from its formation location. 
Thus planet-disk interactions play a crucial role in the lOPF scenario. Here we present numerical 
simulations that first assess the degree of migration for planets of various masses that are forming at 
the DZIB of an active accretion disk, where the effective viscosity is undergoing a rapid increase in 
the radially inward direction. We find that torques exerted on the planet by the disk tend to trap 
the planet at a location very close to the initial pressure maximum where it formed. We then study 
gap opening by these planets to assess at what mass a significant gap is created. Finally we present 
a simple model for DZIB retreat due to penetration of X-rays from the star to the disk midplane. 
Overall, these simulations help to quantify both the mass scale of first, “Vulcan,” planet formation 
and the orbital separation to the location of second planet formation. 

Subject headings: protoplanetary disks, planet-disk interactions, planets and satellites: formation 


1. INTRODUCTION 

Since launch in 2009, Kepler has revealed more t han 
4000 exoplanet candidates (e.g., Mullally et al.|20i5 ). A 
large percentage (> 30%) of these candidates are in sys- 
tems with tightly-packed inner planets (STIPs). These 
are systems that usually have 3 or more detected planets 
of radii ~ 1 — 10 an d with periods less than 100 days 


(Fang & Margot 2012). There are two main scenarios 
that can produce such close-in planets: (1) form ation in 


NelsonI 2012 Gossou et al. 

2013 

20141): (2) formation in 

situ (Hansen & Murray 

20 

m 20TH 

Chiang & Laughlin 

2013 Chatterjee & 'i’an 

20 

14| hereafter 


ets that are trapped in orbits of low order mean motion 
resonances, which is not a particular feat ure of STIPs 
( Baruteauet al.|20T4 Fabrycky et al.|20T4 ). Thus it has 


been proposed that the lack of resonant pile-ups might 
be explained by a lower efficiency of resonance trapping 
or breaking of resonance by later dynamical proc esses for 


the typically low-mass Kepler-detected planets (Goldre- 

ich fc Schlichting||2014 Chatterjee fc Ford||2015L 

'i'he in situ formation scenario faces the chalL 


le nge of 

concentrating enough solids in the inn er region (|Ray- 


mond & Cossou 2014[ [Schlichting 2014). For example, 
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supply of pebbles by radial drift may be truncated if 
planet formatio n, perhaps initiated by the pe bble stream- 
ing instability (lYoudin fc Goodman||2005|) , occurs in 
the outer disk (iLamb rechts V Johansen|2(J14 Levison et 


al. 2015 Bitsch et al. 2015 b In situ formation models 
also face the challenge ot reproducing the observed mass 
versus orbital radius distributions once effects of gas on 
protoplanet migrat ion during the oligar chic growth phase 
are accounted for (|Ogihara et al.||2015 1. 

The Inside-Out Hlanet Formation (lUPF) scenario pro¬ 
posed by CT14 is a new type of in situ formation model. 
It starts with pebble delivery to the midplane transi¬ 
tion region between the innermost MRI-active zone and 
a nonactive “dead zone,” where there is a local pressure 
maximum. The pebbles trapped by the pressure maxi¬ 
mum will build up in a ring, which then forms a proto¬ 
planet, perhap s involving a variety of pr ocesses includ- 
ing streaming ( Youdin fc Goodman||2005 ), gravitational 
( Toom re||1964|) and/or Ross by wave (IVarniere & Tagger] 
2lJU6|) instabilities. The protoplanet is expected to con¬ 


tinue its growth, especially by pebble accretion, until it 
becomes massive enough to open a gap in the disk. This 
gap pushes the pressure maximum outwards by a few 
Hill radii thus creating a new pebble trap that is dis - 
placed from the planet’s orbit (Lambrechts et alj2014) 

"be 


but may also allow MRI-activation in the region beyond 
the planet, which could induce further outward retreat 
of the DZIB. 


As discussed by CT14, planetary migration, either be¬ 
fore gap opening (Type I) or after gap opening (Type 
H), may in principle alter this scenario. However, the ex¬ 
pectation is that the pressure maximum associated with 
this inner disk transition region will also act as a planet 
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trap (Masset et al.||2006 

Matsumura et al.||2009 

Kretke 

& Lin 


Bitsch et al. 

2014), so the planet will stay 


materials, at least until the point of gap opening. How¬ 
ever, this needs to be confirmed for the particular disk 
structure involved in lOPF and this is one of the goals 
of this paper. 


The disk studied by Masset et al. (2006 1 was a purely 


passive disk, i.e, the tempera ture was set only by irra¬ 
diation from the central star. Zhang et al. (2014) stud¬ 
ied planetary migration in a disk with an inner region 
heated by viscous dissipation and an outer region heated 
by stellar luminosity. This structure also contained a 
local pressure maximum that acted as a trap for Type 
I planetary migration. In the model for lOPF, the Hrst 
planet is formed at the outer edge of the MRI active zone 
around the protostar (where T ^ 1200 K, leading to ther¬ 
mal ionization of alkali metals) and migration needs to 
be studied in the context of an active disk, i.e., heated by 
viscous accretion, with accretion rates ~ 10 “® Mq yr“^. 

The magnitude of the planet mass that leads to gap 
opening and potential DZIB retreat is also of crucial im¬ 
portance f or the lOPF theory. Gap op ening has been 


studied by Lin & Papaloizou (1986 
“viscous-thermal” criterion tha' 


1993), who proposed 


sets a gap opening 


mass. This mass is sensitive to the local disk viscosity. 
Magnetohydrodynamic simulations with realistic turbu¬ 
lence have also be en carried out to study gap opening 


in tur bulent disks (|Nelson & Papaloizou|2003 


Zhu et al. 


The second goal of this paper is to study gap open¬ 
ing and quantify the gap opening mass in the context 
of the lOPF disk model, i.e., an inner region of an ac¬ 
tive disk where the effective viscosity is undergoing a 
sharp increase in the radially inward direction. These 
gap opening masses can then be compared to masses 
of the innermost, so-called “Vulcan” planets. An ini¬ 
tial comparison of the simple, analytic gap opening mass 
predicti on with the observ e d Vul cans has been carried 
out by Chatterjee & Tan (2015 hereafter CT15, Pa¬ 
per II), who found a predicted scaling and normaliza¬ 
tion of planet mass versus orbital radius, r, of Mp = 
Mg ~ 5 . 0 (/)g,o. 3 (’'/ 0.1 AU)a _3 Mq, where 4>g, 0.3 is a di¬ 
mensionless parameter that indicates the fraction of the 
Lin-Papaloizou mass scale that is assumed to lead to a 
deep enough gap to truncate planetary accretion of peb¬ 
bles and initiate DZIB retreat. The normalization of this 
gap opening planet mass at given radius also depends on 
the Shakura-Sunyaev viscosity parameter, a, in the dead 
zone inner boundary region, with the fiducial value of 
10“^ being adopted in the above formula. However, this 
value is quite uncertain and the observed planet masses 
may require a somewhat lower value of 2 x 10“^ (Paper 
H), or a lower value of (pc- 

Finally, our third goal is investigate how gap opening 
may induce DZIB retreat, which then sets the location 
of second planet formation. Here we will present a sim¬ 
ple, heuristic first exploration of this process by modeling 
the location of the DZIB as set by penetration of X-rays 
emitted from the protostellar corona to the disk midplane 
beyond the planet-induced gap. This location then sets 
the radius for an imposed radial profile of viscosity that 
simulates the transition region from an MRI-active inner 
region to an outer dead zone. Increased viscosity leads to 


reduced densities that make it easier for X-rays to pen¬ 
etrate further. We present example toy models of this 
process in which the location of the DZIB ends up sta¬ 
bilizing several initial gap widths away from the planet. 
Such separations can be compared to the observed sepa¬ 
ration of orbits of innermost STIPs planets. 

In we describe an analytic model to calculate our 
disk parameters. In we describe our numerical set-up 
and test the numericm simulations against the analytic 
model. In (Q we study migration of planets of various 
masses located in the DZIB region. In ^ we study gap 
opening and DZIB retreat. In ^ we discuss the impli¬ 
cations of our results for observed planets. We conclude 
in 


2. ANALYTIC ESTIMATES OF DISK STRUCTURE AND 
GAP OPENING MASS 


As in Paper I, we follow Frank et al. (2002 )’s deriva¬ 
tion of the structure of a viscously heated accretion disk. 
These results will be compared to the numerical simu¬ 
lations with FARGO, described below in To achieve 
consistency with these simulations we fina we need to 
make a small correction in the choice of the vertical op¬ 
tical depth equation and now adopt: 


T = 0.5Sk, 


( 1 ) 


with the factor of 0.5 being introduced since the disk 
has two faces (or equivalently a change in the definition 
of “midplane” conditions). Thus when calculating the 
balance between energy dissipation and radiative cooling, 
the optical depth is integrated from the midplane to each 
surface of the disk. This change then leads to modest (< 
40%) changes in the normalization of the disk structure 
equations compared to CT14. For example, for disk mass 
surface density we now derive: 
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(CTI4 derived a normalization of 106 gem ^), where 
/r = 2.33?71h = 3.90 x 10 “^^ g is the mean particle 
mass (assuming nue = O.lnu), fcs is Boltzmann’s con¬ 
stant, 7 = 1 . 471,4 is the power law exponent of the 
barotropic equation of state P = Kp^ where we have 
normalized for H 2 with rotational modes excited, osb is 
Stefan-Boltzmann’s constant, m* = is the stel¬ 

lar mass, K = KiqIO cm^ g“^ is disk opacity (normalized 
to ex pected protoplanetary disk values, e.g.,|Wood et al.| 


20021 , /r = 1 — ^JrPJr, (where r* is stellar radius), and 


m = TO-glO ® Mq yr ^ is the accretion rate. 
For the disk aspect ratio we find 


h 

r 


X a-'/'" (/.m)’''" C'” (3) 


The above equations will be used to set the initial con¬ 
ditions of the FARGO simulations described below. 
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Given this disk structure, the viscous criterion for the 
gap opening planet mass is : 

40^771* 


Mg = 4>g - 


I'K 


(4) 


Implementing our disk model, we obtain: 


Mq = 40(/)g 


128/ 
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-2/5 ( M 


(f) 


-4/5 
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4/5 -1/5 


^SB 


X a4/5G-7/10^3/10^l/5 (J^^)2/5 ^1/10 ( 5 ) 

^ 14.08(/)G7f.'4 


This is a factor 0.745 smaller than the CT14 result. 

The full set of disk structure equations with our revised 
normalizations are presented in Appendix [X} 


3. NUMERICAL SIMULATION METHODOLOGY 


We used the 2D code FARGO-ADSG jBaruteau & 


|2008ail: 
Tout with 


Massetl 2008ajb I, which is built on FARGu (iVlasset 


(p 


2 UUUp, but with the energy equation and disk self-gravity 
Implemented. Self-gravity has not been turned on in our 
calculation since the disk is far from gravitational insta¬ 
bility in our problem, which focusses on inner regions 
near the star. We also implemented a simple method of 
radiative cooling that uses the disk surface density to cal- 
culate the optica l depth and cooling rate as described in 
Zhu et al. (2012 1 . We set up a disk with the same mass 
surface density profile, flaring index, opacity, adiabatic 
index, mean molecular mass and central stellar mass as 
the fiducial analytic model in ^ 

The energy equation imple mented in FARGO-ADSG 
is ( Baruteau fc Masse^|2008a ): 


dc -* 

— -I- V • {ev) = - (7 - 1) eV • F -I- (5+ - Q_, ( 6 ) 


where e is the thermal energy per unit area, v is the flow 
velocity and Q+{Q-) denote heating (cooling) source 
terms, assumed to be positive quantities. In a steady 
accretion disk, these terms can be written as 


SgT ks 

7 - 1 F 


Gm*- -Tin 

7-1 





dvi 

dxj 


= VTg 



Q- 


IScTSB 

3t 


^-1 


(7) 

( 8 ) 

(9) 


where S'i j are components of the viscous stress tensor, 

T is the midplane temperature. Note t = O.bTgK is 
the optical depth from disk midplane to the surfac^ 
First we set up an accretion disk with a constant value 
of a, i.e., without a transition zone, and a steady ac¬ 
cretion rate of 10 “® Mq yr“^, simulating a range of 
radii from 0.02 to 0.3 AU. In all simulations we use 
the EVANESGENT boundary condition. This damps 
the disk values (surface density, velocities and energy 
density) to the initial axisymmetric disk conditions. 
The damping regions are rings within the radial ranges 


5 During numerical testing we found an error in the public ver¬ 
sion of the code: it only calculates half the value of vYigidv^jdr)'^. 


[fmin: l-25rniin] and [0.84rmaxj Linax]; where 7‘iaiin(fmax) IS 
the inner (outer) radius of the disk. The damping pa¬ 
rameter is increased from zero to a maximum from the 
edge of damping zone to the edge of the disk. 

Within each orbit, the damping amplitude is not a 
constant: the actual damping is the calculated damp¬ 
ing function times a coefficient, which gradually grows 
from zero to one, linearly with time. 

For our standard “medium resolution (MR),” the disk 
is evenly divided by 300 sectors azimuthally and loga¬ 
rithmically divided by 128 sectors radially. This gives a 
typical grid size of 0.002 x 0.002 AU, which is about the 
same as the Hill radius of a l-Mc-mass planet located at 
0.1 AU. We also carry out “high resolution (HR)” simu¬ 
lations that have twice the MR resolution, and a few test 
runs (of more limited time duration) at “super-high res¬ 
olution” (SHR) that have four times the MR resolution. 

We next set up a disk that has a radial jump in a. 
Moving inwards, a rises from odzib = 0.001 at 0.1 AU 
to aMRi = 0.01 at 0.07 AU, i.e., over a transition width 
of At-dzib = 0.03 AU, with the transition described by 
part of the function odzib + (1 + sinx)(aMRi — Q^dzib) 
from X = 37 r /2 (at 0.1 AU) to a: = tt/2 (at 0.07 AU). The 
initial mass surface density profile is also chosen to give 
a steady accretion rate across this transition region. The 
steady radial profile achieved after 1000 orbits is close to 
this initial choice and is shown in Fig. 



r(AU) 

Figure 1. Mass surface density profile after 1000 orbits of our 
fiducial “transition” disk where a rises from 0.001 at 0.1 AU to 
0.01 at 0.07 AU with steady accretion rate of 10“® Mq yr“^. The 
dashed line is the profile of a “constant o” disk (i.e., a = 0.001). 

4. PROTOPLANET MIGRATION 

We study the migration of protoplanets of various 
fixed masses (0.1, 0.5 and 1.0 Mq) that are inserted 
into the transition zone region of the disk by measur¬ 
ing the torques exerted on the planet by the disk. For 
each planet mass, we run simulations holding the planet 
at a constant radius, exploring uniformly from 0.085 to 
0.115 AU with a spacing of 0.001 AU. So for each planet 
mass there are 31 MR and HR simulations run, for a total 
of 186 simulations. Note we keep the a viscosity radial 
profile constant in all the simulations in this section. We 
run each simulation for 400 orbits, thus allowing the gas 
disk to achieve a quasi equilibrium structure. We then 
evaluate the torque on the planet by averaging over the 
next 100 orbits. Given the 2D nature of these simula¬ 
tions, we adopt a smoothing length in the calculation of 
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Figure 2. Transition zone migration torque, F, profiles for 
0.1 Mq (top), 0.5 Mq (middle) and 1.0 Mq (bottom) planets 
held at fixed radius. Medium resolution (MR) and high resolution 
(HR) runs are shown, with the torque induced by the disk on the 
planet being the average from orbit 400 to 500 after introduction 
of the planet. A positive torque leads to outward migration; a 
negative torque to inward migration. The HR results show F de¬ 
creasing from positive to negative values in the transition region 
(where a starts decreasing as r increases just interior of 0.1 AU) 
. We expect the planet will migrate inwards from positions where 
the total torque is negative and migrate outwards where the total 
torque is positive. It will stop at zero torque. Thus these torque 
profiles are indicative of there being a stable “ planet trap ” in this 
location, where F = 0. 


the p lanet’s gravitational potential of 0.6h (Miiller et al 


20121 . 


The torque, T, is scaled as Tq = {Mp/mt.y'Egr‘^flph~^, 
where r is the orbital radius of the planet and Eg is the 
average gas mass surface density at this orbital radius. 
The grid size in our standard MR run is more than twice 
the Hill radius of a planet with Mp = O.lMc- We there¬ 
fore carry out HR simulations with twice greater resolu¬ 
tion and examine numerical convergence. 

In all cases for 0.1 Me, 0.5 Mq and 1.0 Mo-mass 
planets (which correspond to very different levels of gap 
opening—see below), the torque profiles show a common 


behavior of T decreasing from positive to negative values 
in the transition region (where a starts decreasing as r 
increases just interior of 0.1 AU), indicative of a stable 
“planet trap” in this location, where T = 0 (see Figurej^. 
We note that the magnitude of T decreases significanHy 
as Mp increases from the Type I to Type H regimes. 

This confirmation of planet trapping at the DZIB tran¬ 
sition zone is a key requirement of lOPF, since this al¬ 
lows a planet to continue to grow by pebble accretion 
from low to relatively high masses. Next we investigate 
how pebble accretion may be disrupted by gap opening 
by looking at the detailed structure of the gaps. 

5. GAP OPENING 

5.1. Gap Structure & Opening Mass for Fixed a Profile 

The gap opening process in lOPF is critical to termi¬ 
nation of pebble accretion of the first planet, thus setting 
its mass. It may also be important, in combination with 
dead zone retreat, in setting the location of the pressure 
maximum that leads to new pebble ring formation and 
then second planet formation. 

We now investigate the disk structure that is induced 
by introducing planets of various masses at a fixed loca¬ 
tion of r = 0.1 AU, i.e., at the transition zone of the disk 
with a fixed a profile. Figure shows the radial profiles 
of mass surface density, midplane pressure and midplane 
temperature after 1000 orbits of evolution from intro¬ 
duction of planets of masses Mp = 0.1 Mq to 1.0 Mq in 
steps of 0.1 Ma- We note that the gas profile settles very 
quickly, within ^ 100 orbits, to a profile close to that of 
the final state at 1000 orbits. 

As planet mass increases, the radial profiles show a 
gradual deepening of the gap in both the mass sur¬ 
face density and pressure profiles. The temperature also 
dips due to reduced viscous heating and smaller opti¬ 
cal depths. For this particular set-up, with a transition 
zone width of 0.03 AU, it is only once planet masses 
are > O.bMo that there is significant displacement of 
the azimuthally averaged pressure maximum away from 
the planet’s orbital radius. Table lists the separa¬ 
tion of the pressure maximum from the planet’s or¬ 
bital radius at 0.1 AU, Arp^naxi but normalized by the 
planet’s Hill radius, Rh = (Mp/{3mt))^^^r. We define 
0Ar,Pmax = ^Tpniax/RH- As Mp increases from 0.4 to 
0.5 Ma, (/'Ar.Pmax grows by a factor of 10 from ~ 0.4 to 
~ 5. Figure shows this also corresponds to the mass 
surface density maximum retreating outwards by a sim¬ 
ilar amount. 

We have compared gap opening results for simulations 
with Mp = 0.1,0.5, l.OMfj at MR and HR resolutions. 
We find very similar radial profiles of E and P, with 
agreement in values at better than 10% across the width 
of the gap region. 

Figure]^ shows the 2D view of the perturbation of the 
disk mass surface density (compared to the disk with no 
planet) in the r-cji plane. The strengthening spiral arms, 
along with deepening gap, are evident. 

5.2. DZIB Retreat with an Evolving a Profile 

When a planet opens a gap in the disk leading to de¬ 
creasing density in the vicinity of the planet, we expect 
higher ionization levels given the dependence of re¬ 
combination rates, and thus greater likelihood of acti¬ 
vating the MRI, which then results in the retreat of the 
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Table 1 

Retreat of the local pressure maximum during gap opening 


MpIMq 

1 0.1 

0.2 

0.3 

0.4 

0.5 

0.6 

0.7 

0.8 

0.9 

1.0 

fliAr.Pmax = Rh 

1 0.61 

0.48 

0.42 

0.38 

4.75 

4.47 

5.09 

4.86 

4.68 

5.27 


^The separation between new pressure maximum and planet orbital radius, scaled to planet’s Hill radius. These results are for disks with 
fixed Q profile with a transition width of ArDZiB = 0-03 AU. Note, location of pressure maximum is measured at the center of grid cell of 
maximum pressure. 



0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 
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Figure 3. Radial structure of disk mass surface density (top), 
midplane pressure (middle) and midplane temperature (bottom) 
during gap opening for planets of mass Mp = 0.1 Mq to 1.0 Mg in 
steps of 0.1 Me held fixed at the transition zone radius of 0.1 AU. 
Profiles are shown for medium resolution runs after 1000 orbits. 
Gap opening manifests itself via a decrease in S, P and T in the 
vicinity of the planet and, once Mp > 0.5Mq, a significant re¬ 
treat of the local surface density and pressure maxima outwards 
by several Hill radii. 


dead zone inner boundary. Higher levels of ionization 
may also arise due to shock heating in spiral arms that 
are induced by the planet. We note that these spiral 
arms, propagating outwards into the dead zone, will by 
themselves also provide a source of enhanced effective 
viscosity. The lower densities of the gap region may also 
allow for increased penetration of X-rays that raise the 
ionization level and activate the MRI. Such effects could 
cause the viscosity to make its radial transition further 
out in the disk, which would lower the density in the gap 
region and vicinity even more. 

The full treatment of these processes is a very complex 
problem, beyond the scope of the present paper. Here we 
present a simple, heuristic model for DZIB retreat based 
on X-ray penetration from the protostellar corona to the 
disk midplane. The X-rays are assumed to propagate 
from a height of 0.05 AU at r = 0 (i.e.. Hi?©, which is 
expected to be ~ 3-4 stellar radii for a solar mass pre- 
main sequence star). The thickness of the transition zone 
from the MRI-active region to dead zone is assumed to be 
characterized by i-e., a mass surface density that 

scales with the X-ray penetration mass surface density 
due to absorption. Ex- 

Starting from the FARGO simulation with (j)c = 0.5 
after 300 orbits (i.e., when a wide, but relatively shallow, 
gap has just been opened), the radius of the DZIB, rDziB, 
is re-evaluated each orbit by the condition that the mass 
surface density along the path from the X-ray corona to 
the midplane equals (px^x- For simplicity, we assume 


there is negligible material along this path length at lo¬ 
cations interior to the planet, i.e., at r < 0.1 AU, which 
we consider to be a good approximation due to both vis¬ 
cous clearing of the inner disk and also since the disk’s 
vertical scale height is small in this region. Thus our 
path integral starts from r = 0.1 AU and extends out to 
■cdzib- Note, that this path starts with a relatively shal¬ 
low angle to the midplane of ~ 24°, which then decreases 
as the disk evolves and tozib increases. 

In addition, we need to specify the structure of the 
a transition zone in the X-ray penetration region. We 
expect it to be broader than the transition zone set by 
thermal ionization of alkali metals. For simplicity, we 
assume the functional shape of the radial variation of 
a is the same as that adopted in The outer radius 
where a starts rising from the DZIB value of 0.001 is set 
by the mass surface density along the pathlength being 
equal to (^x^x- We note that the width of this tran¬ 
sition region may be somewhat larger than the actual 
penetration depth of the X-rays (i.e., > 1), because 

of, for example, the outward propagation o f turbulence, 
which has b een seen in the simulations of |Dzyurkevich| 


et al. (2010). For the inner radius, r^, where a reaches 
the tiili MrI value of 0.01, we investigate three cases, i.e., 
inner radii of 0.07 AU (Case A: the same as adopted for 
our migration and gap opening studies, above), 0.1 AU 
(Case B; the location of planet), and 0.11 AU (Case C: 
the location of the initial gap outer edge that has been 
induced by the planet). 

For a simple estimate on the values of Ex and (f>x 
to be used in this model, we note that X-ray opacity, 
even at ^ 10 keV, may be dominated by heavy ele¬ 
ments in the gas, rather than in dust, especially if dust 
has mostly coagulated into large grains and pebbles that 
have settled to the midplane and been trapped at the 
DZIB pressure maximum. For interstellar gas with so¬ 
lar a bundances and th e case of heavy elements deple¬ 
tion, Igea & Glassgold (1999) find an absorption opac¬ 
ity of te = 0.056(E/1 g cm“^)(i?/10 keV)“^-®^ (valid 
for E = 2-30 keV). Thus for nokeV = 1,10 requires 
Ex = 18,180 gcm“^ and raokeV = 1,10 requires Ex = 
39, 390 gcm“^. Thus we choose a characteristic value 
of Ex = lOOgcm”^, equivalent to absorption optical 
depths of a few for ~ 20 keV X-rays. We note, however, 
that the radiative transfer of these higher energy X-rays 
is likely to be controlled by Compton scattering by elec¬ 
trons present in neutral H, H 2 and He, which may allow 
increased penetration to the disk midplane at greater ra¬ 
dial distances via scattering from the disk atmosphere. 
For this reason, along with the extra thickness of the 
transition zone that is expected due to propagation of 
turbulent wakes from the MRI-active zone (including ver¬ 
tically from modest scale heights), we choose a fiducial 
value of 4>x = 5. 

From the above discussion, it is clear that more de- 
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Figure 4. 2D r-<j) view of the mass surface density perturbation, E/Sq, induced by different planet masses from 0.1 to 1.0 Mq (where 
So is the surface density in the absence of a planet). 



Figure 5. Radius of DZIB (i.e., the outer transition radius of a 
profile) as set by X-ray penetration in evolving a. disks. Medium 
resolution (MR) results are shown with dashed lines; high resolu¬ 
tion (HR) results with solid lines. The three cases of inner tran¬ 
sition radii of 0.07 AU (Case A, red lines), 0.1 AU (Case B, blue 
lines) and 0.11 AU (Case C, green lines) are shown. In each case, 
the estimated converged, final retreat radius is shown with a hori¬ 
zontal dot-dashed line (see text). 


tailed calculations are needed to better constrain the 
combined parameter (j)x'^x^ including ioni zation equilib- 
rium calcul ations, e.g., similar to those of|Igea fc Glass-] 


gold (1999) and then applying the resulting resistivities 
to determine the extent of the MRI-active zone (Mohanty 
& Tan, in prep.). The time variable nature of the X-ray 
luminosity may also need to be accounted for. 

Figurej^shows results from these DZIB retreat models. 
We find that increases gradually over timescales 

of several thousand orbits of the planet. For Case A, by 
6500 orbits the rate of increase is very slow, i.e., from the 
6500th to 7000th orbit the percentage increase is <0.06% 
for the MR run and 0.07% for the HR run. The equiva¬ 


lent results for Cases B and C are similar. However, the 
value of rDziB after 7000 orbits does depend somewhat on 
the resolution of the simulation: for the Case A MR run 
?’DZiB = 0.128 AU; for the HR run cdzib = 0.139 AU, i.e., 
the amount of retreat (Arpmax) is about 38.6% greater 
in the HR run. We thus also carried out twice higher, 
SHR, resolution runs for ^1000 orbits. We find that 
Arpniax is ^5% greater than in the HR simulation. We 
use these results to estimate that the final numerically 
converged retreat distance is ~5-10% greater than the 
HR result (or about 40% greater than the MR result), 
i.e., rDziB(CaseA) ~ 0.145 AU, i.e., Arpmax — 0.045 AU. 
For this case with Mp = 0.5 Mg = 5.59 M®, for which 
the planet’s Hill radius is Rh = 1.78 x 10“^ AU then 
<?i’Ar,Pmax = Arp„iax/7?_f/ = 25.3, which is a significantly 
greater retreat (i.e., ^ 5x greater) than due to simple 
gap opening in the fixed a model. 

For Cases B and C the amount of retreat is larger and 
takes somewhat longer to stabilize (although the retreat 
has essentially stopped by 7000 orbits). After 7000 or¬ 
bits, the HR simulations find a 24.0% and 20.8% larger 
retreat than the MR runs for Cases B and C, respectively. 
The SHR results at ~ 1000 orbits are 5.6% and 4.9% 
greater, respectively. Thus we estimate final converged 
values are about 25% greater than the MR results, i.e., 
for Case B DZIB retreat to 0.20 AU (Arpmax — 0.10 AU; 
<('Ar,Pmax = 56.3) and for Case C DZIB retreat to 0.22 AU 
(Arpinax — 0.12 AU; (/lAr.Pmax = 67.6). These results are 
also summarized in Tabled 

As a check on self-consistency, we have also investi¬ 
gated the potential effect of X-ray penetration inducing 
an evolving a profile in a disk that does not yet have a 
planet or gap. For this model, where the accretion disk 
extends uninterrupted to very close to the protostar at 
~ 0.02 AU, we set the height of the stellar X-ray corona 
to be equal to this inner disk radius, i.e., 0.02 AU. We 
also start the path integral for X-ray penetration from 
this distance. We find that for the same choices of Ex 
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and (j)x adopted above, the X-rays do not penetrate be¬ 
yond ^0.1 AU and so the DZIB a profile set by thermal 
ionization of alkali metals is not expected to be affected. 
However, we find that utilizing a larger X-ray corona 
height of 0.05 AU does enable X-ray penetration to be¬ 
yond 0.1 AU in the no-planet disk. Thus, in the context 
of this simple model, DZIB retreat is also influenced by 
properties of the scaleheight of X-ray emission from the 
protostar. 

Overall, we therefore see that the amount of DZIB re¬ 
treat, which in the lOPF scenario is expected to set the 
location of second planet formation, depends somewhat 
sensitively on the detailed structure of the a transition 
zone and the X-ray emission properties of the protostar. 
Thus future work to model this zone including the full 
physics of activation of the MRI in response to a changing 
ionization fraction, e.g., due to X-rays or other ionization 
sources, is needed. However, we can still conclude that, 
for the parameters explored here, the amount of retreat 
is significantly greater than the retreat due to initial gap 
opening. In ^we compare these retreat distances to the 
orbital spacings of the Kepler-ohseived planets. 

Figure shows the radial structure of the disks for 
these evolving a models after 7000 orbits, and compares 
to the fixed a profile model, which was already settled 
after 1000 orbits. This figure shows that the local pres¬ 
sure maximum remains closely located with the outer a 
transition radius, i.e., where it begins to rise from its 
DZIB value of O.OOI. 

Figure shows the 2D r-0 plots of perturbations of 
mass surface density and pressure for the Case A evolv¬ 
ing a model and compares to the fixed a model, both 
of which have Mp — 0.5 Mq. Figure shows the abso¬ 
lute values of these quantities. For completeness, to bet¬ 
ter illustrate physical appearance of disk structures, we 
also show true spatial maps of disk mass surface density 
structure of Case A in Figure!^ These figures reveal the 
decreasing absolute mass surface density of the induced 
spiral arm in the case where DZIB retreat has proceeded 
further. 

We emphasize that we have presented a very simple, 
heuristic model for a disk in which the viscosity profile 
responds actively to X-ray penetration (which itself is 
sensitive to the disk’s structure). More accurate mod¬ 
eling would involve self-consistent calculation of the vis- 
cosity caused by the MRI in the transi tion zone (e.g., 
Bai & Stone||201I Mohanty et al. 2013), including cal- 
culation of the ionization structure of the disk as it is 
affected by gap opening. The results that are needed 
from such modeling are the shape and absolute values of 
the disk viscosity in the MRI transition zone, which could 
then be implemented parametrically in, e.g., FARGO hy¬ 
drodynamic simulations of the disk. However, potential 
instability of the dead zone i nner edge, which may lead 
it to vary its radial location (Latter &: Balbusl|20 I2 ), or 
vortex formation at the inner edge (I'aure et al. 2DI5), 
should also be assessed. Although the model presented 
so far should be regarded as illustrative, as we discuss in 
the scale of pressure maximum retreat that it exhibits 
is comparable to the orbital separations of the innermost 
planets in STIPs. 


6. IMPLICATIONS FOR OBSERVED PLANETS 


CTI4 examined the separation of planet orbits in 
STIPs (normalized by Hill radius of inner planet in the 
pair), finding a broad range of separations from inner¬ 
most, Vulcan, planet to next planet of 4>Ar,i "^few-IOO, 
with a broad peak from ~ 20 — 60. Note, however, that 
this analysis was based on the simple, single power law 
planetary mass-size relation of Lissauer et al. (2011). 

CTI5 analyzed the available mass-size data tor STIPs 
planets and derived an improved, piece-wise power law 
mass-size relation. Here, we use this relation (PL3 in 
CTI5) to re-evaluate the distributions of (jiAr,! through 
4 >Ar,4y where the numerical value of the subscript refers 
to the number of the planet that is the innermost in the 
pair being considered, countin g ou t from the star. These 
results are presented in Figure [T0| 

Note some of the dispersion in separations is likely to 
be induced by an intrinsic dispersion in the densities of 
the planets, leading to inaccurate estimation of masses. 
For example, the observed dispersion in density of a fac¬ 
tor of about 5 (e.g., CTI5) leads to a dispersion in in¬ 
ferred mass of the same factor and a dispersion in Rh 
and cj)Ar of a factor of 1.7. Also, the observed values of 
4 >Ar,i rnay be overestimated if some fraction of second 
planets are non-transiting and thus not detected due to 
slight orbita l m isalignments. 

In Figure 10 we also show the locations of (/>Ar.Pmax 
for the fixed and evolving a disk models we considered. 
We note that the evolving a disk models can achieve 
separations of ~ 20 to 70 i?//, which overlap with the 
observed distribution of 

As discussed by CTI4, the distributions of ()>Ar,i are 
significantly different from those of the outer separa¬ 
tions. This conclusion remains unchanged by our use 
of the improved piecewise power-law CTI5 mass-size re¬ 
lation. For KPC systems with Np > 3, the probability 
that ((>Ar,i and (/>Ar ,2 are drawn from the same distribu¬ 
tion is 9 X 10“^. Similarly, for systems with Np > 4, 
the respective probabilities that two (pAr distributions 
are drawn from the same underlying distribution are 
7 X I0“®, 2 X I0“®, and 0.36 for ^Ar.i and 4 ’ Ar , 2 , '('Ar.i 
and 4 >Ar,3, and 4 >Ar,2 and 4 >Ar,3- This indicates that the 
distribution of separations between the first two planets 
is different than between the 2nd and 3rd and other outer 
pairs of planets in these systems. 

As also discussed by CTI4, the generally larger values 
of (j)Ar,i may be a signature of the greater relative effect 
of inner disk clearing after first, Vulcan planet formation 
on the location of the second planet, compared to later, 
more incremental dead zone retreats. 


7. SUMMARY AND DISCUSSION 

We have explored three main aspects of the Inside-Out 
Planet Formation (lOPF) scenario with numerical sim¬ 
ulations of planet-disk interactions. Our main findings 
are the following: 

(1) Planets of all relevant masses are trapped at the 
dead zone inner boundary transition region, so planetary 
migration (of Type I or Type H) does not appear to be a 
problem for lOPF. Especially during the early stages of 
potential Type I migration, the fact that the protoplanet 
is trapped at its formation location should allow it to 
continue to accrete, especially by pebble accretion. 

(2) The mass scale of significant gap opening that can 
affect pebble accretion by displacing the local pressure 


















Table 2 

DZIB retreat for evolving a models for different inner transition radii, ri 


Models 

Medium Res. (7000 orbits) 

^Ar,Pmax 

High Res. (7000 orbits) 

Estimated Convergence 

Case A (0.5Mq, ri=0.07 AU) 
Case B (0.5Mq, rj=0.10 AU) 
Case C (0.5M(5, ri=0.11 AU) 

15.9 (0.028 AU) 

41.2 (0.073 AU) 

52.4 (0.093 AU) 

22.1 (0.039 AU) 

51.1 (0.091 AU) 

63.3 (0.112 AU) 

~25.3 (0.045 AU) 
~56.3 (0.100 AU) 
~67.6 (0.120 AU) 


^ ^Ar,Pmax measures the separation between the planet (at 0.1 AU) and the outer transition of a profile (as set by X-ray penetration) in 
units of the planet’s Hill radius. The outer transition radius is at a location very similar to that of the local pressure maximum (see Fig.[^. 
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Figure 6. Radial structure (Left panel: mass surface density; Right panel: midplane pressure) of a disk with a 0.5 Mq planet held fixed 
at 0.1 AU. 1st row: The fixed a profile disk model (Q:(r) is shown by the dotted red line). Disk structure is shown after 1000 orbits (solid 
blue line), compared to initial state (dashed green line). The extent of the planet’s Hill sphere is shown by the vertical shaded band. 2nd 
row: The evolving a disk model (HR run) at 7000 orbits due to X-ray penetration in which the inner a transition radius is held fixed at 
0.07 AU (Case A). 3rd row: Same as above, but for Case B with inner transition radius at 0.1 AU. Bottom row: Same as above, but for 
Case C with inner transition radius at 0.11 AU. 


maximum away from the planet is Mp ~ O.^Mq^ i.e., 
(j)G — 0.5, which is about 6 Mq in the fiducial case at 
r = 0.1 AU. This is similar to the fiducial value adopted 
by CT14, especially since given the renormalization of 
Mg by a factor of 0.745 in this paper (§2), this means 
that the overall change in planet mass is only a factor of 
1.24. However, these results are derived in the context 
of pure hydrodynamic simulations and may be sensitive 
to including extra physics of MRI activation in the gap 
region. 

(3) A simple model for MRI activation in the gap re¬ 
gion, implemented by an evolving a viscosity profile set 
by X-ray penetration, leads to greater retreat of the pres¬ 
sure maximum out to separations of ^ 20 to 70 Hill radii 
of the planet, but sensitive to model assumptions about 
the shape of the viscosity profile in the MRI-active region 
to dead zone transition region. Such separations overlap 
with the observed orbital separations of innermost plan¬ 


ets in STIPs. Future work is needed to better model 
viscosity profiles in such transition regions. 

While the above results are supportive of the lOPF 
model for being relevant for formation of STIPs, a num¬ 
ber of open questions remain. For example, even though 
planets with masses > O.IMg appear likely to be trapped 
at a location where they can continue to grow by pebble 
accretion, it is possible that this is n ot the case at earlier 
stages for lower mass protoplanets. Faure et al. (2015) 
have pointed out that inwardly migrating vortices can 
form at the dead zone inner edge and they can poten¬ 
tially interact with the planet, possibly causing it to also 
migrate inwards. However, their 3-D simulation results 
are for an unstratified disk and the mass accretion rate 
is not a constant across the dead zone inner edge, lead¬ 
ing to mass pile up and potentially influencing vortex 
formation. The potential effects of vortex interactions 
still needs to be confirmed in 3-D stratified simulations 
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Figure 7. Mass surface density (top row) and midplane pressure (bottom row) perturbation plots of disks with a 0.5 Mq planet with 
(Left column) fixed a profile (after 1000 orbits) and (Right column) evolving a profile (Case A, after 7000 orbits). 


where the presence of the surface active layer may trans¬ 
port the disk mass smoothly across the dead zone inner 


planetesimals combined with continued pebble accretion. 
Improved study of pebble supply truncation during the 


the dead zone region to represent the effects of an active 

tions (see Lambrechts et al. 

2014 

layer and/or turbulent wakes spreading from the inner 

tions focussed on the outer c 

usk), e 


2014 for example calcula- 


MRI-active region, which may control the accretion rate 
in this radial region of the disk. With this set-up we 
do not observe the sharp density pe aks at the dead zone 


inner edge that were present in the Fame et al. (20151 
simulations. 

Further work is needed on the lOPF model to study the 
transition of the pebble ring to a single dominant proto¬ 
planet, potentially involving streaming and gravitational 
instabilities of the pebble population and/or Rossby wave 
instabilities in the gas, followed by oligarchic growth of 


Potential migration of planets after their main accre¬ 
tion phase and retreat of their natal DZIB, including due 
to interactions between neighboring planets, is another 
topic for future study. However, we note that for the 
fiducial lOPF model, the gas mass that is present in the 
vicinity of first planet to form is relatively small com¬ 
pared to the planet’s mass, thus limiting the scope of its 
migration. 


























10 



0.05 0.10 0.15 0.20 0.05 0.10 0.15 0.20 

r (AU) r (AU) 


Figure 8. Same as Figure but now showing absolute values of mass surface density (top row) and midplane pressure (bottom row) of 
disks with a 0.5 Mq planet \rah {Left column) fixed a profile (after 1000 orbits) and {Right column) evolving a profile (Case A, after 7000 
orbits). 
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viding computational resources and support that have 
contributed to the research results reported in this pub¬ 
lication. 
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Figure 9. Same as Figure[^top row, but now showing a spatial map of the disk structure induced by a 0.5 Mq planet with (Left) fixed 
a profile (after 1000 orbits) and (Right) evolving a profile (Case A, after 7000 orbits). 



Figure 10. A comparison of orbital separations of the Kepler-detected planets (see legend) and our simulation results for DZIB pressure 
maximum retreat: vertical long-dashed line shows (j5>Ar,Pmax = 4.75 for the initial gap-openning in the fixed a profile disk. The three short 
dashed lines show the estimated values of 0Ar,Pmax = 25.3,56.3,67.6 for Cases A, B, C, respectively, of the evolving a profile disk. 
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APPENDIX 

DERIVATION OF DISK PROPERTIES 


Here we update the analysis presented in CT14, which involves modest changes in disk properties, as discussed in 
^ Conservation of angular momentum in a viscous disk implies: 


dt 

where r(r, t) is the shear torque from viscosity: 




27r dr 


T{r,t) = 2TTrvY,„r‘^ —. 

dr 


In a steady disk, we have ^ = 0, so integrating 


Al 


we obtain: 


Substituting T(r,t), we obtain: 


r'Enr^flVr = - -h C. 

27r 


dfl C _ 

,.3 “ 


The boundary condition at r = r* + 5r, ^ = 0, so constant C is: 


C = + 5r) = —m(G'm*r*)^/^/(27r), 


where 


m = 2TTrTig{—Vr)- 

For Sr r^, we have 12 (r* + Sr) « 12;^ (r*) = Gm*/ r^, and substituting into 


A4 


we have: 


m 

= to 




1 / 2 - 


So the energy dissipation rate per unit area (only via one face of the disk) at radius r is: 


D{r) = 


r{dil,/dr) fdn\^ iGm^fn 

Anr 2 \ dr J Sirr^ 


'-(7) 


1/2 


(Al) 

(A2) 

(A3) 

(A4) 

(A5) 

(A6) 

(A7) 

(AS) 
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From the mid-plane to disk surface, the optical depth is: 

not = O.SEgK 


(A9) 


From Eq. 3.11 in Hubeny (1990), we know the relation between surface temperature and actual temperature where 
the optical depth is r to the surface of the disk: 


T^ir) = -X^T 1 - 


2ni 


— -T^ T 

o eff • 


Here Ttot = 0.5Ek, and in our situation, D{r) = Tc = T (not), so 


16cr 

3Sk 


= Dir). 


We have an expression for mid-plane temperature: 


4 _ 9Gmi,m'SgK 
° 1287r err 3 


'-(7) 


1/2' 


but we still have an unknown E profile, so now we try to deduce it. In an a disk, we have: 

^2 


V = OL- 


UJk 


Substituting in A7 we obtain: 


_ m^/Gml 
® 37rac^ 


1 - 


\ 1/2' 


For an adiabatic sound speed Cg, we have: 


Then the mid-plane pressure satisfies: 


(7) 

2 dP P 

Cg = -^ = 7- 

dp p 


p ^ pfcTe _ pkTc 


,.-3/2 


p 3c 

Now we have the following relation between Eg, r and Tc'. 

ihy/Gm^ p 




® 37ra "fkTc 


'-(7) 


1 / 2 ' 


,,-3/2 


so, substituting for Tc{r) from equation A12 the surface density profile is: 

and mid-plane temperature (here we use T instead of Tc) is: 


(AlO) 

(All) 

(A12) 

(A13) 

(A14) 

(A15) 

(A16) 

(A17) 

(A18) 

(A19) 


For scale height h, usually it is assumed that the disk is vertically isothermal in a passive disk, here in an active disk, 
we use an adiabatic vertical structure: 

h cg Cg cg ^^20) 


r VK yjGm^lr 


so the aspect ratio is: 




(A21) 


Since the disk vertical structure satisfies p{z) = p{z = 0)e mid-plane density p = Tigl^/^h can be expressed 


as: 




(A22) 




























